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Abstract 

Motivated by finite element spaces used for representation of temperature in the compatible fi¬ 
nite element approach for numerical weather prediction, we introduce locally bounded transport 
schemes for (partially-)continuous finite element spaces. The underlying high-order transport 
scheme is constructed by injecting the partially-continuous field into an embedding discontinuous 
finite element space, applying a stable upwind discontinuous Galerkin (DG) scheme, and project¬ 
ing back into the partially-continuous space; we call this an embedded DG transport scheme. We 
prove that this scheme is stable in provided that the underlying upwind DG scheme is. We 
then provide a framework for applying limiters for embedded DG transport schemes. Standard 
DG limiters are applied during the underlying DG scheme. We introduce a new localised form of 
element-based flux-correction which we apply to limiting the projection back into the partially- 
continuous space, so that the whole transport scheme is bounded. We provide details in the specific 
case of tensor-product finite element spaces on wedge elements that are discontinuous Pl/Ql in 
the horizontal and continuous P2 in the vertical. The framework is illustrated with numerical 
tests. 

Keywords: Discontinuous Galerkin, slope limiters, flux corrected transport, 
convection-dominated transport, numerical weather prediction 


1. Introduction 


Recently there has been a lot of activity in the development of finite element methods for 
numerical weather prediction (NWP), using continuous (mainly spectral) finite elements as well as 


discontinuous finite elements ([Fournier et ah, 

2(1(141 

Thomas and Loft 

2005t 

Dennis et ah 

2011 

Kelly and Giraldo, 2012 

[ Giraldo et ah, 2013 [Marras et ah, 2013t Brdar et ah, 20131 Bao et ah 

2015 

); see 

Marras et ah ( 

2015 

) for a comprehensive review. A key aspect of NWP models is the neec 


for transport schemes that preserve discrete analogues of properties of the transport equation such 
as monotonicity (shape preservation) and positivity; these properties are particularly important 
when treating tracers such as moisture. Discontinuous Galerkin methods can be interpreted as 
a generalisation of finite volume methods and hence the roadmap for the development of shape 


preserving and positivity preserving methods is relatively clear (see Gockburn and Shu (2001) for 


an introduction to this topic). However, this is not the case for continuous Galerkin methods, 
and so different approaches must be used. In the NWP community, limiters for GG methods have 


been considered by Marras et ah (2012), who used first-order subcells to reduce the method to 


first-order upwind in oscillatory regions, and Guba et ah (2014), who exploited the monotonicity of 


the element-averaged scheme in the spectral element method to build a quasi-monotone limiter. 
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In this paper, we address the problem of finding snitable limiters for the partially continuous 
finite element spaces for tracers that arise in the framework of compatible finite element methods 


for numerical weather prediction models ( 

Cotter and Shipton 

2012 

Gotter and Thuburn, 

2014 

Staniforth et al. 

2013 

McRae and Gotter 

2014 

). Gompatible finite element methods have been 


proposed as an evolution of the C-grid staggered finite difference methods that are very popular 
in NWP. Within the UK dynamical core “Gung-Ho” project, this evolution is being driven by the 
need to move away from the latitude-longitude grids which are currently used in NWP models, 
since they prohibit parallel scalability (Staniforth and Thuburn, 2012). Compatible finite element 
methods rely on choosing compatible finite element spaces for the various prognostic fields (velocity, 
density, temperature, etc.), in order to avoid spurious numerical wave propagation that pollutes the 
numerical solution on long time scales. In particular, in three dimensional models, this calls for the 
velocity space to be a div-conforming space such as Raviart-Thomas, and the density space is the 
corresponding discontinuous space. Many current operational forecasting models, such as the Met 


Office Unified Model (Davies et ah, 2005), use a Charney-Philips grid staggering in the vertical, 
to avoid a spurious mode in the vertical. When translated into the framework of compatible 
finite element spaces, this requires the temperature space to be a tensor product of discontinuous 
functions in the horizontal and continuous functions in the vertical (more details are given below). 
Physics/dynamics coupling then requires that other tracers (moisture, chemical species etc.) also 
use the same finite element space as temperature. 

A critical requirement for numerical weather prediction models is that the transport schemes 
for advected tracers do not lead to the creation of new local maxima and minima, since their 
coupling back into the dynamics is very sensitive. In the compatible finite element framework, this 
calls for the development of limiters for partially-continuous finite element spaces. Since there is 


a well-developed framework for limiters for discontinuous Galerkin methods (Biswas et al. 


Burbeau et ah, |2001 

Gockburn and Shu, 12001; 

Hoteit et al. 

and Aliabadi 

2005 

I 

Kuzmin 

2010 

Zhang and Shu 

2011 

), in 


2004 ; Krivodonova et al. 2004t Tu 


approach of (i) injecting the solution into an embedding discontinuous finite element space at 
the beginning of the timestep, then (ii) applying a standard discontinuous Galerkin timestepping 
scheme, before finally (iii) projecting the solution back into the partially continuous space. If the 
discontinuous Galerkin scheme is combined with a slope limiter, the only step where overshoots 
and undershoots can occur is in the final projection. In this paper we describe a localised limiter 


for the projection stage, which is a modification of element-based limiters (Lohner et ah, 1987 


Kuzmin and Turek, 2002) previously applied to remapping in Lohner (2008); Kuzmin et al. (2010). 


This leads to a locally bounded advection scheme when combined with the other steps. 
The main results of this paper are: 


1. The introduction of an embedded discontinuous Galerkin scheme which is demonstrated to 
be linearly stable, 

2. The introduction of localised element-based limiters to remove spurious oscillations when 
projecting from from discontinuous to continuous finite element spaces, which are necessary 
to make the whole transport scheme bounded, 

3. When combined with standard limiters for the discontinuous Galerkin stage, the overall 
scheme remains locally bounded, addressing the previously unsolved problem of how to limit 
partially continuous finite element spaces that arise in the compatible finite element frame¬ 
work. 


Our bounded transport scheme can also be used for continuous finite element methods, although 
other approaches are available that do not involve intermediate use of discontinuous Galerkin 
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methods. 

The rest of the paper is structured as follows. 


particular, more detail on the hnite element spaces is provided in Section 2.1 


The problem is formulated in Section In 

The embedded 


discontinuous Galerkin schemes are introduced in Section 2.2 it is also shown that these schemes 
are stable if the underlying discontinuous Galerkin scheme is stable. The limiters are described in 


Section 2.3 In Sectionwe provide some numerical examples. Finally, in Section]^ we provide a 
summary and outlook. 


2. Formulation 

2.1. Finite element spaces 

We begin by defining the partially continuons finite element spaces nnder consideration. In 
three dimensions, the element domain is constructed as the tensor prodnct of a two-dimensional 
horizontal element domain (a triangle or a qnadrilateral) and a one-dimensional vertical element 
domain (i.e., an interval); we obtain triangular prism or hexahedral element domains aligned with 
the vertical direction. For a vertical slice geometry in two dimensions (freqnently used in testcases 
dnring model development), the horizontal domain is also an interval, and we obtain quadrilateral 
elements aligned with the vertical direction. 

To motivate the problem of transport schemes for a partially continnous finite element space, 
we first consider a compatible finite element scheme that uses a discontinuous finite element space 
for density. This is typically formed as the tensor prodnct of the DGk space in the horizontal 
(degree k polynomials on triangles or bi-fc polynomials on qnadrilaterals, allowing discontinnities 
between elements) and the DGi space in the vertical. We consider the case where the same degree 
is chosen in horizontal and vertical, i.e. k = I, althongh there are no restrictions in the framework. 
We will denote this space as DGk x DGk- 

In the compatible finite element framework, the vertical velocity space is staggered in the 
vertical from the pressnre space; the staggering is selected by reqniring that the divergence (i.e., 
the vertical derivative of the vertical velocity) maps from the vertical velocity space to the pressnre 
space. This means that vertical velocity is stored as a field in DGk x GGk+i (where GGk+i denotes 
degree k + 1 polynomials in each interval element, with G^ continuity between elements). To avoid 
spnrious hydrostatic pressnre modes, one may then choose to store (potential) temperatnre in the 
same space as vertical velocity (this is the finite element version of the Gharney-Phillips staggering). 
Figure [T] provides diagrams showing the nodes for these spaces in the case k = 1. Details of how 
to antomate the construction of these hnite element spaces within a code generation framework 


are provided in McRae et ah (2015) 


Monotonic transport schemes for temperatnre are often required, particularly in challenging 
testcases snch as baroclinic front generation. Fnrther, dynamics-physics conpling reqnires that 
other tracers snch as moistnre must be stored at the same points as temperature; many of these 
tracers are involved in parameterisation calculations that involve switches and monotonic advection 
is reqnired to avoid spnrious formation of rain patterns at the gridscale, for example. Hence, we 
must address the challenge of monotonic advection in the partially continuous DGk x GGk+i space. 

In this paper, we shall concentrate on the case of DGi x GG 2 . This is motivated by the fact 
that we wish to use standard DG npwind schemes where the advected tracer is simply evalnated 
on the upwind side; the lowest order space DGo x GGi leads to a first order scheme in this case. 
We may retnrn to higher order spaces in fnture work. 
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Figure 1: Diagrams showing the nodes for the various spaces in the compatible finite element framework in the 2D 
vertical slice case and degree k = 1. From left to right: the velocity space RTi, the pressure space DGi x DGi, the 
vertical part of the velocity space RTi, and the temperature space DGi x CG 2 . Circles denote scalar nodes, whilst 
lines denote normal components of a vector. 
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Figure 2: Diagrams showing the nodes for the partially continuous space V and the discontinuous space V^ in the 
case V = DGi x GG 2 and V = DGi x DG 2 . 
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2.2. Embedded Discontinuous Galerkin schemes 

In this section we describe the basic embedded transport scheme as a linear transport scheme 
without limiters. The scheme, which can be applied to continuous or partially-continuous finite 
element spaces, is motivated by the fact that limiters are most easily applied to fully discontinuous 
hnite element spaces. We call the continuous or partially-continuous finite element space V, and 
let V be the smallest fully discontinuous finite element space containing V. A diagram illustrating 
V and V in our case of interest, namely the hnite element space for temperature described in the 
previous section, is shown in Figure 

Before describing the transport scheme, we make a few dehnitions. 

Definition 1 (Injection operator). For u & V <Z V, we denote I \ V ^ V the natural injection 
operator. 

The injection operator does nothing mathematically except to identify lu as a member of V 
instead of just V. However, in a computer implementation, it requires us to expand m in a new 
basis. This can be cheaply evaluated element-by-element. 

Definition 2 (Propagation operator). Let A : V V denote the operator representing the ap¬ 
plication of one timestep of an L'^-stable discontinuous Galerkin discretisation of the transport 
equation. 

For example, A could be the combination of an upwind discontinuous Galerkin method with a 
suitable Runge-Kutta scheme. 

Definition 3 (Projection operator). For u & V we define the projection P : V ^ V by 

{v,Pu) = {v,u), Vn G V. 


In a computer implementation, this requires the inversion of the mass matrix associated with 

R. 

We now combine these operators to construct our embedded discontinuous Galerkin scheme. 

Definition 4 (Embedded discontinuous Galerkin scheme). Let V <ZV, with injection operator I, 
projection operator P and propagation operator A. Then one step of the embedded discontinuous 
Galerkin scheme is defined as 

gn+l ^ pAjgn^ QU^ QU+l ^ 


The L^ stability of this scheme is ensured by the following result. 

Proposition 5. Let a > 0 be the stability constant of the the propagation operator A, such that 

\\Az 


= sup 
zev,\\z\\>o 


< a. 


( 1 ) 


where || ■ || denotes the norm. Then, the stability constant a* of the embedded discontinuous 
Galerkin scheme on V satisfies a* < a. 


Proof. 


\\PAIz\\ \\PAIz\\ 

sup —^— = sup II . ,| 
^ey,||2||>o IGII 2ey,||2||>o 


< sup 
2ey,||2||>o 


\\PAz\\ 


< sup 
2£y,||2||>o 


||Hz| 


< a. 


( 2 ) 


as required. In the last inequality we used the fact that \\Pz\\ < ||;§||, which is a consequence of 
the Riesz representation theorem. □ 
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Corollary 6 . For a given velocity field u, let A^t denote the propagation operator for timestep 
size At. Let At* denote the critical timestep for A^t, i-e., 

II^Atll < 1, for At < At*. 

Then, the critical timestep size AA for the embedded discontinuous Galerkin scheme PA^th is at 
least as large as At*. 

Proof. If At < At*, then 

WpAaJW < II^Atll < 1, 

as required. □ 


Hence, the embedded DG scheme is stable whenever the propagation operator A is. 

For the numerical examples in this paper, we consider the case V = DGi x CG 2 (our temper¬ 
ature space) and V = DGi x DG 2 . For a given divergence-free velocity held u, dehned on the 
domain G and satisfying w ■ n = 0 on the domain boundary dfl, A represents the application of 
one timestep applied to the transport equation 


e^ = -u-ve = -V -{uO), 


(3) 


discretised using the usual Runge-Kutta discontinuous Galerkin discretisation (see Gockburn and 
for a review). To do this, hrst we dehne L : R —)■ R by 


Shu 


( 2001 ) 


'yL6dx = —At / ■ u6 dx + At / [it 7 ] 6 *dS', 


(4) 


JO JO Jr 

where F is the set of interior facets in the hnite element mesh, with the two sides of each facet 
arbitrarily labelled by -|- and —, the jump operator denotes ■ n+ + v~ ■ n~, and where 9 

is the upwind value of 6 dehned by 


0 


0 + if w ■ n+ < 0 , 
9~ otherwise. 


Then, the timestepping method is dehned by the usual 3rd order 
method (Shu and Osher, 1988), 


step oorrtiv Limesieppmg 


= 9^ + L9^, 

(5) 

02 ^ _gn _^0l 2^01^^ 

( 6 ) 

AQn ^ gn +1 ^ ^ 

3 3 

(7) 


Since the hnite element space V is discontinuous in the horizontal, the projection P : R —>■ R 
decouples into independent problems to solve in each column [i.e., the mass matrix for DGi x GG 2 
is column-block diagonal). 


2.3. Bounded transport 

Next we wish to add limiters to the scheme. This is done in two stages. First, a slope limiter 
should be incorporated into the R propagator. A; we call the resulting scheme A. A suitable 

After replacing A with A, the only way that the solution can 


2.3.1 


limiter is dehned in Section 
generate overshoots and undershoots is after the application of the projection P. To control these 
unwanted oscillations, we apply a (conservative) hux correction to the projection, referred to as 
hux corrected remapping (Kuzmin et ah, 2010); this is described in Section 2.3.2 We denote 


the hux corrected remapping P, and the resulting bounded transport scheme may be written as 

gn+l ^ 
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2.3.1. Slope limiter for the propagator A 

In principle, any suitable discontinuous Galerkin slope limiter can be used in the propagator A. 
In this paper we used the vertex-based slope limiter of Kuzmin (2010). This limiter is both very 
easy to implement, and supports a treatment of the quadratic structure in the vertical. Before 
presenting the limiter for V = DGi x DG 2 (recall that this is the space we must use to obtain a 
transport scheme for our DGi x GG 2 space used for temperature), we first review the concepts in 
the simpler case oiV = DGi x DGi. The basic idea for 6 & V = DGi x DGi is to write 


e = e + Ae, 


( 8 ) 


where 6 is the projection of 6 into DGq, i.e. in each element 6 is the element-averaged value of 
6. Then, for each vertex i in the mesh, we compute maximum and minimum bounds and 

by Computing the maximum and minimum values of 6 over all the elements that contain that 
vertex, respectively. In each element e we then compute a constant 0 < c^e < 1 such that the value 
of 

^mm,i A: T ^ ^max,i) (9) 


at each vertex i contained by element e. The optimal value of the correction factor Oe can be 
determined using the formula of Barth and Jespersen (1989) 


Op = min 
ieA4 


min 11 , s' 
1 


min 


1 , 



if > 0 , 

if = 0 , 

if < 0 , 


( 10 ) 


where Me is the set of vertices of element e and is the unconstrained value of 

the DGi shape function at the i-th vertex. 

For our temperature space DGi x DG 2 applied to numerical weather prediction applications, 
we assume that we have a columnar mesh. This means that the prismatic elements are stacked 
vertically in layers, with vertical sidewalls (but possibly with tilted top and bottom faces to facili¬ 
tate terrain-following meshes, so that the elements are trapezia). This allows us to adopt a Taylor 
basis in the vertical, i.e. the basis in local coordinates is the tensor product of a Taylor basis in 
the vertical with a Lagrange basis in the horizontal. We write 


e = e + {ei-e) + {e-ei), 


( 11 ) 


where 6*1 G DGi x and satisfies the following conditions: 

1 . h = 0, 

2 . ^ and take the same values along the horizontal element midline in local coordinates. 

Then, § e DGi x DGi whilst ^ G DGi x DGq. 

First, we limit the quadratic component in the vertical (the third term in Equation ([IT|)), 
performing the following steps. 

1. In each element, compute and evaluate the derivative at the horizontal cell midline to 
obtain ^ G DGi X DGq. If the quadratic component 6 — 61 is limited to zero then ^ will 
become equal to 

2. In each column, at each vertex i, compute bounds |min,i and |max,i by taking the maximum 
value of ^ at that vertex in the elements sharing that vertex in the column. 
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3. In each element, compute element correction factors ai^e according to 


min<'l, 


Ml _Ml 

^zlmax,^ gzle,i 


_M| 

dze,i = * 


ai e = mm < 

’ ZGjVe 


1 

min ^ 1, 


Ml _Mi 

|min,i 

'I .-S»| - 

,le,z 32 |c,* 


if 

if 

if 


d® I 
d® I 

g^le,z 

mi 

dzle,i 


Ml 

dz\ 

Ml 

dz\ 

Ml 

dz\ 


e,i > 0, 

e,i d; 
e < 0. 


(12) 


This approach can also be extended to meshes in spherical geometry for which all side walls are 
parallel to the radial direction^ having replaced by the radial derivative. 

Second, we apply the vertex-based limiter to the DGi x DGi component ^i, obtaining limiting 
constants uq. We then finally evaluate 


9 1-^ 9 = 9 + ao{9i — 9) + ai{9 — 9i). 


(13) 


To reduce diffusion of smooth extrema, it was recommended in Kuzmin (2010) to recompute the 
ao coefficients according to 

Q;o,e H- max(Q!o,e, ai,e)- (14) 


However, this does not work in the case of DGi x DG 2 since there is no quadratic component 
in the horizontal direction, and hence nonsmooth extrema in the horizontal direction will not be 
detected. A possible remedy is to use a;o,e for the horizontal gradient and max(a!o,e 5 for the 
vertical gradient or to limit the directional derivatives separately using an anisotropic version of 
the vertex-based slope limiter (Kuzmin et ah (2015)). 

This limiter is applied to the input to A and after each SSPRK stage, to ensure that no new 
maxima or minima appear in the solution over the timestep. 


(Kuzmin et ah (2015 


2.3.2. Flux corrected remapping 

The final step of the embedded DG scheme is the projection P of the DG solution (which we 
denote here as 9) back into V. We obtain a high-order, but oscillatory solution, which we denote 
9^. To obtain a bounded solution, we introduce a localised element-based limiter that blends 9^ 
with a low-order bounded solution 9^, such that high-order approximation is preserved wherever 
overshoots and undershoots are not present. 

First, we must obtain the low-order bounded solution. Using the Taylor basis, we remove the 
quadratic part of 9, to obtain 9 G DGi x DGi. A low-order bounded solution can then be obtained 
by applying a lumped mass projection, 

P m 

Mi9^ = / (j)i9dx = ^Qik9k, i = l,...,m, (15) 

k=i 

where the lumped mass M is defined by 

Mi= f (pidx, (16) 

Jq 

the projection matrix Q is defined by 

Qik= / (j)ii/Jkdx, i = 1,... ,n,/c = 1,... ,m, (17) 

Jn 


^Such meshes arise when terrain following grids are used in spherical geometry. 














1 is a Lagrange basis for DGi x DGi and {4>i}'^^i is a Lagrange basis for DGi x GG 2 - 
The Inmped mass M and projection matrix Q both have strictly positive entries. This means 
that for each 1 < i < n, the basis coefficient Of is a weighted average of valnes of 9 coming from 
elements that he in S'(z), the snpport of 0^. The weights are all positive, and hence the valne of Of 
is bonnded by the maximnm and minimnm valnes of 9 in S{i). Hence, no new maxima or minima 
appear in the low order solntion. 

Next, we combine the low order and high order solntions element-by-element, in a process 


called element-based flux correction. Element based flux correction was introduced in (Lohner 


et al., 1987) and formulated for conservative remapping in (Lohner, 2008 Kuzmin et ah, 2010). 


Here, we use a new localised element-based formulation, where element contributions to the low 
and high order solutions are blended locally and then assembled. 


To formulate the element-based limiter, we note that the consistent mass counterpart of (15) 
is given by 


'^Mij9f = (j)i9dx, i = l,...,n, 

Jn 


i=i 


where 


Mij — 


(l)i(j)j dx. 


(18) 


(19) 


'n 


First, by repeated addition and subtraction of terms, we write (with no implied sum over the 
index i) 

M,9^ = M,9^ + U 


where 


h = M,9f - Mij9f + M,9f + 


H 

j ’ 


= M,9, 


H 


W MijO^ + / 4>i{0 - 9) dx. 
, Jn 


( 20 ) 

( 21 ) 

( 22 ) 


This can be decomposed into elements to obtain 

MiO? = E (M'Of + ft ). ft = M'Of - Y. MJOf + juo- 0) d X, 

e J 

where 

Mf = y^jdx, and = JcpiCpjdx. 

Importantly, the contributions /f of element e to its vertices sum to zero, since 

n n n n p n 

Yf‘=T, ^‘‘0^ -EE + / E 

i=l i=l i=l 7=1 i=l 


(23) 


(24) 


= 1 


n n p 

E ‘''‘to" - E + / (« - 9) d 

*=1 i=i _ 


X = d. 


(25) 


=0 


=0 
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It follows that the total mass of the solution remains unchanged (i.e., ) 

if all contributions of the same element are reduced by the same amount. 

We can then choose element limiting constants a® to get 

= Y. • P6) 

e 

where 0 < Oe < 1 is a limiting constant for each element which is chosen to satisfy vertex bounds 
obtained from the nodal values of 6. 

The bounds in each vertex are obtained as follows. First element bounds 9^^^^ and 9 ^^„ are 
obtained from 9 by maximising/minimising over the vertices of element e. Then for each vertex i, 
maxima/minima are obtained by maximising/minimising over the elements containing the vertex: 


9max,i 

e e 

The correction factor Oe is chosen so as to enforce the local inequality constraints 


< AffOf + aj‘ < M‘e, 


(27) 


(28) 


Summing over all elements, one obtains the corresponding global estimate 

(29) 

which proves that the corrected value := 9^ + ^ Yhe bounded by 9 m^^ i and 9 m^r, i- 

To enforce the above maximum principles, we limit the element contributions /f using 


Mi9min,i < Mi9f + 


aef! < MA 


max,I ) 


CKe = min 
i&Me 


min <1 
1 


L ’ fr 

min j 


l ’ f! 


if /f>0, 
if /f = 0, 
if /f<0. 


(30) 


This definition of cxe corresponds to a localised version of the element-based multidimensional FCT 
limiter ((Lohner et ah, 1987 Kuzmin and Turek, 2002)) and has the same structure as formula (10) 
for the correction factors that we used to constrain the DGi approximation. A further advantage 
of the localised formulation is that the limited fluxes can be built independently in each element, 
before assembling globally and dividing by the global lumped mass by iterating over nodes. 


3. Numerical Experiments 

In this section, we provide some numerical experiments demonstrating the localised limiter for 
embedded Discontinuous Galerkin schemes. 


3.1. Solid body rotation 

In this standard test case, the transport equations are solved in the unit square D = (0,1)^ 
with velocity field u{x, y) = (0.5 — y,x — 0.5), i.e. a solid body rotation in anticlockwise direction 
about the centre of the domain, so that the exact solution at time f = 27r is equal to the initial 
condition. The initial condition is chosen to be the standard hump-cone-slotted cylinder configu¬ 
ration defined in LeVeque (1996|, and solved on a regular mesh with element width h = 1/100 and 
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Figure 3: Solid body rotation test (Section [3d] ) on a lOOx 100 grid. Solution is interpolated to a DGi discontinuous 
field before plotting. Left: initial condition. Right: Solution after one rotation. The solution is free from over- 
and undershoots and exhibits comparable numerical dissipation to discontinuous Galerkin methods combined with 
limiters. 


Courant number 0.3. The result, shown in Figure]^ is comparable with the result for the DGi dis¬ 
continuous Galerkin vertex-based limiter shown in Figure 2 of Kuzmin (2010); it is free from over- 


and undershoots and exhibits a similar amount of numerical diffusion. It is also hard to distinguish 
between the a:-direction, where the finite element space is discontinuous, and the t/-direction, where 
the space is continuous. This suggests that we have achieved our goal of constructing a limited 
transport scheme for our partially-continuous finite element space. 

3.2. Advection of a discontinous function with curvature 

In this test case, the transport equations are solved in the unit square G = (0,1)^ with velocity 
held u = (1,0), i.e. steady translation in the x-direction (which is the direction of discontinuity 
in the hnite element space). The initial condition is 


9 = 


4|/(1 — y) + I if 0.2 < X < 0.4, 
4|/(1 — y) otherwise. 


(31) 


This test case is challenging because the height of the “plateau” next to the continuity varies as a 
function of y {i.e., in the direction tangential to the discontinuity); this means that the behaviour 
of the limiter is more sensitive to the process of obtaining local bounds. 

The equations are integrated until t = 0.4 in a 100 x 100 square grid and Courant number 0.3. 
The results are showing in Figure]^ One can see qualitatively that the degradation in the solution 
due to the limiter and numerical errors is not too great. 


3.3. Convergence test: deformational flow 

In this test, we consider the advection of a smooth function by a deformational flow field that 
is reversed so that the function at time t = 1 is equal to the initial condition. As is standard for 
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Figure 4: Results from the test case in Section |3.2| Left: The initial condition. Right: The solution at time 
t = 0.4. 


this type of test, we add a translational component to the flow and solve the problem with periodic 
boundary conditions to eliminate the possibility of fortuitous error cancellation due to the time 
reversal. 

The transport equations are solved in a unit square, with periodic boundary conditions in the 
x-direction. The initial condition is 

9{x, 0) = 0.25(1 + cos(r)), r = min (^0.2, ^y{x - 0.3^ + {y - 0.5^/0.2^ , 
and the velocity field is 

u{x, t) = {1 — 5(0.5 — t) sin(27r(a; — t)) cos(7r?/), 5(0.5 — t) cos(27r(x — t)) sin(7rj/)), 

where x = {x,y). The problem was solved on a sequence of regular meshes with square elements 
at hxed timestep At = 0.000856898, and the error was computed. A plot of the errors is 
provided in Figure As expected, we obtain second-order convergence (the quadratic space in 
the vertical does not enhance convergence rate because the full two-dimensional quadratic space 
is not spanned). 


4. Summary and Outlook 


In this paper we described a limited transport scheme for partially-continuous finite element 
spaces. Motivated by numerical weather prediction applications, where the finite element space for 
temperature and other tracers is imposed by hydrostatic balance and wave propagation properties, 
we focussed particularly on the case of tensor-product elements that are continuous in the vertical 
direction but discontinuous in the horizontal. However, the entire methodology applies to standard 
Co hnite element spaces. The transport scheme was demonstrated in terms of convergence rate on 
smooth solutions and dissipative behaviour for non-smooth solutions in some standard testcases. 

Having a bounded transport scheme for tracers is a strong requirement for numerical weather 
prediction algorithms; the development of our scheme advances the practical usage of the com¬ 
patible finite element methods described in the introduction. The performance of this transport 
scheme applied to temperature in a fully coupled atmosphere model will be evaluated in 2D and 
3D testcases as part of the “Gung Ho” UK Dynamical Core project in collaboration with the Met 
Office. In the case of triangular prism elements we anticipate that it may be necessary to modify 


the algorithm above to limit the time derivatives as described in Kuzmin (2013). 
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Ax 

error 

0.05 

0.02 

0.0125 

0.01 

0.0319911 

0.0048104 

0.0017125 

0.0010108 


Figure 5: Convergence plot for deformational flow experiment (Section |3.3[ ) showing second order convergence, and 
a table of error values. 
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A key novel aspect of our transport scheme is the localised element-based FCT limiter. This 

limiter has much broader potential for use in FCT schemes for continuous hnite element spaces, 

which will be explored and developed in future work. 
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